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Abstract 

In a neutron lifetime experiment conducted at the National Institute of Standards 
and Technology, protons produced by neutron decay events are confined in a proton 
trap. In each run of the experiment, there is a trapping stage of duration r. After 
the trapping stage, protons are purged from the trap. A proton detector provides 
incomplete information because it goes dead after detecting the first of any purged 
protons. Further, there is a dead time 5 between the end of the trapping stage in 
one run and the beginning of the next trapping stage in the next run. Based on the 
fraction of runs where a proton is detected, I estimate the trapping rate A by the 
method of maximum likelihood. I show that the expected value of the maximum 
likelihood estimate is infinite. To obtain a maximum likelihood estimate with a fi- 
nite expected value and a well-defined and finite variance, I restrict attention to a 
subsample of all realizations of the data. This subsample excludes an exceedingly 
rare realization that yields an infinite- valued estimate of A. I present asymptotically 
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valid formulas for the bias, root-mean-square prediction error, and standard devia- 
tion of the maximum likelihood estimate of A for this subsample. Based on nominal 
values of A and the dead time 5, I determine the optimal duration of the trapping 
stage r by minimizing the root-mean-square prediction error of the estimate. 

Key words: , Lifetimes, 21. 10. Tg, Nuclear tests of fundamental interactions and 
symmetries, 24. 80+y, Properties of protons and neutrons, 14.20. Dh, Probability 
theory, stochastic processes, and statistics,02.50-r 



1 Introduction 

Ion traps play a key role in fundamental physics experiments (Ref. 1). In this 
paper, I focus on statistical methods for uncertainty analysis and planning of 
proton trap neutron lifetime experiments (Refs. 1-5) and related experiments 
such as Ref. 6. When a neutron decays, it produces a proton, an electron and 
an antineutrino. An accurate determination of the mean lifetime of the neutron 
is critically important for testing the fundamental theories of physics (Ref. 7). 
Further, the mean lifetime of the neutron is an important parameter in the 
astrophysical theory of big bang nucleosynthesis (Ref. 8). In a proton trap neu- 
tron lifetime experiment performed at the National Institute of Standards and 
Technology (NIST), a beam of neutrons passes through a detection volume. 
Based on measurements of the neutron flux and the proton production rate, 
one measures the mean lifetime of the neutron. Each run of the experiment 
consists of trapping stage where protons are confined in a trap (Refs. 2-5), 
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and a detection stage. The detector provides incomplete information because 
it goes dead after detecting the first proton. Based on the number of runs 
where a proton is detected, one can estimate the proton trapping rate. 

In earlier work [Ref. 9], this estimation problem was studied using a Bayesian 
method. Given a particular realization of the data (the number of runs where 
at least one ion (proton in this paper) is trapped), a formula for the posterior 
mean of the ion trapping rate was presented based on a prior probability 
model for the trapping rate. In this work, I estimate the trapping rate by the 
method of maximum likelihood and focus on the statistical properties of this 
estimate. I neglect physical sources of systematic error due to effects such as a 
time varying proton trapping rate or fluctuations in the actual trapping stage 
interval about the nominal value sought by the experimenter. 

In Section 2, 1 demonstrate that the bias (expected value minus true value) and 
variance of the maximum likelihood estimate of the trapping rate A are infinite. 
This is so because a rare realization of the data yields an infinite estimate 
of A. This technical problem can be dealt with in various ways. One could 
quantify uncertainty by constructing confidence intervals of finite width even 
though the variance of the estimate is infinite. Another approach would be to 
introduce a stopping rule so that the experiment is continued until no protons 
are trapped in at least one run. I do not pursue either of these approaches here. 
Instead, I restrict the sample space to include only realizations of data where 
one observes at least one run where no protons are trapped. For realizations 
of data in this subsample, the maximum likelihood estimate has finite first 
and second moments. In Section 3, I derive asymptotically valid formulas for 
the bias, variance, and mean-square-error of a maximum likelihood estimate 
of the proton trapping rate computed from this subsample. In general, one 
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expects estimates that are nonlinear functions of the observed data, such as 
the maximum likelihood estimate of the trapping rate, to be biased [Ref. 10]. 
In Section 4, where, based on nominal values of trapping rate and dead time, I 
determine the trapping time that minimizes the root-mean-square prediction 
error of the maximum likelihood estimate of A in the subsample of interest. 



2 Statistical Model 



In a simulated proton trapping experiment there are many runs. During each 
run, I assume that the duration of the proton trapping stage r is an adjustable 
constant that is known with negligible uncertainty. During the trapping stage, 
I assume that protons are trapped at a constant rate A. Further, I restrict 
attention to the case where A > 0. After the trapping stage, protons are purged 
from the trap. A proton detector provides incomplete information because it 
goes dead after detecting the first of any purged protons. Further, there is a 
fixed dead time 5 between the end of the trapping stage in one run and the 
beginning of the next trapping stage in the next run. I assume that 5 is known 
with negligible uncertainty. If the total time of the experiment is T, the total 
number of runs is 

N run = INT{ ). (1) 

T + 

Above, the function INT(x) rounds the continuous variable x down to the 
nearest integer. Let n + be the observed number of runs where at least one 
proton is trapped. I model the number of protons trapped during any run as a 
realization of a Poisson process with expected value Ar. Hence, the probability 
that no ion is trapped for a given run is 
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p = exp(-Ar). (2) 
The maximum likelihood estimate of po is 

1 "run 

where n + is the number of runs where at least one proton is trapped. Thus, 
the maximum likelihood estimate of A is 



A = -ilnp = -iln(l-^-). (4) 

' ' ly run 



Since n + is a binomial random variable, the probability that n + = k is P(k), 
where 



Hence, the expected value of the maximum likelihood estimate of A is 



1 Nrun L 

E(X) = -- £p(fc)ln(l- — ). (6) 
Similarly, the expected squared value of the estimate is 



1 N run jL 

£(A 2 ) = ^£P(fc)(m(l-— )) 2 . (7) 

' fc = -''run 

For A > 0, P(N run ) = (1 - po)^™" > 0, and both the expected value (first 
moment) and expected squared value (second moment) of A are infinite. The 
variance of A, VAR(X), is not defined because 

VAR{\) = E(X 2 ) - (£(A)) 2 , (8) 
and both terms on the right hand side of Eq. 8 are infinite. 
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To ensure that both E(\) and E(\ 2 ) are finite, I restrict the sample space to 
realizations of the data where n + < N run . From a practical point of view, this 
means that realizations of data where n+ = N run would be ignored. For neu- 
tron lifetime experiments of current interest, the probability that n + = N run 
is negligible provided that r is judiciously chosen. Hence, this subsampling 
restriction does not significantly affect data collection procedures for neutron 
lifetime experiments of current interest. In this subsample, the discrete prob- 
ability density function for allowed realizations of n + — 0, 1, • ■ • , N run — 1 is 
P*(k), where 



l P(yN rUn j 



For this subsample, the first two moments of the maximum likelihood estimate 
are 



i N run — 1 7, 

E(X) = -- E p.(*)Mi-tH, (io) 

' k_Q -''run 



and 



i N run — 1 7 

E(X 2 ) = ~ 2 E P*(fc)(m(l--— )) 2 . (11) 

fc=0 l*run 



Since the first two moments (Eqns. 10 and 11) of A are finite, the variance of 
A is defined and finite. Next, I present analytical formulas to approximate the 
fractional bias, fractional standard deviation and fractional root-mean-square 
prediction error of the estimate computed for this subsample. 
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3 Asymptotic Approximations for Analysis in Subsample 



In the subsample where n + < N run , I derive asymptotically valid approxima- 
tions for the fractional bias (FBI AS), fractional root-mean-square prediction 
error (FRMS), and fraction standard deviation (FSE) of A where 



FBIAS= E{X ~ X) , (12) 
A 



FRMS = v^ A ~ A ) 2 ; (13) 
A 



and 
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FSE = v E ( A E ( X )) = J (frmS) 2 - (FBI AS) 2 . (14) 
A 



To facilitate analysis of A in the subsample I write 



Po = Po ~ e = Po(l - — ), (15) 

Po 



e = po-p . (16) 
Thus, 

\np = ln(p ) + ln(l - w) (17) 
where 

w = —. (18) 
Po 

In the subsample, w takes discrete values in the following interval 
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1- — <w<1-—^ — . (19) 
Po N run p 



The BIAS of the maximum likelihood estimate of A in the subsample is 

BIAS = -T'^Eilnpo) - ln(po)) = -r^^h^l - w)) (20) 
where w is a random variable. 

I derive an asymptotically valid expression for BIAS based on a local approx- 
imation for ln(l — w) in the vicinity of 0. I approximate f(w) = — ln(l — w) 
as a fourth order polynomial f(w) where 



2^4 

7#i in in 

f(w) » /(0) + ^(0) + -/ 2 (0) + -/ 3 (0) + -f\0) (21) 
where f k (0) is the kth derivative of f(w) evaluated at w — 0. Thus, 



Since 



"run 



the central moments /x r = E({n + — E(n + )) r ) are relevant. In the full sample 
where n + < N run , n + is a binomial random variable with an expected value 
equal to N run {l — po). Hence, in full sample, its first four central moments are 

[11] 

fii = 0, (24) 
A*2 = N run p (l -po), (25) 
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A*3 = N run p (l - p )(2po - 1), (26) 

and 

^4 = 3iV r 2 un #j(i _ po )2 + ]V runPo (i _ Po )(i _ 6po (l -p )). (27) 

Since F(rz + = iV run ) tends to exponentially as a function of N run , the asymp- 
totic central moments of n + in the subsample are given by Eqns. 24-27. 

Based on Eqns. 22,25,26 and 27, 1 get the following approximation for FBI AS 



~ \t [ 2(N runPo )^ 3(N runPo )^ 4(N runPo r l ' 1 J 



In a similar calculation where I approximate f(w) as a second order polyno- 
mial, I get the following approximation for FRMS 



~ Ar V (A^nPo) 2 + (iV rMnPo )3 + 4(7V rU nR>) 4 • 1 ) 

From Eqns. 28 and 29, I derive an approximation for FSE 



FSE = \J(FRMS) 2 - (FBI AS) 2 . (30) 



Since the asymptotic standard deviation of the random variable w is 



A-Po 



(31) 



I expect Eqns. 28-30 to be asymptotically valid as N run increases to large 
values. Next, I present evidence consistent with this expectation for a number 
of cases (Tables 1,2). 
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4 Example 



For particular cases, I compute the actual values of FRMS, FBI AS and 
FSE using Eqns. 9-14. I set 5 = 100 fjs (0.0001 s) and A = 1 Hz because these 
are typical values for experiments done at NIST. For experiments of total 
duration of T = 10 s ,25 s, 50 s, 100 s, 200 s, and 400 s, FBI AS was much 
less than FRMS (see Figure 1 and Table 1). Furthermore, FBI AS was more 
sensitive to r than FRMS was (Figure 1). For the cases summarized in Table 
1, the fractional systematic error and fractional RMS prediction error are well 
approximated as FBI AS oc T _1 , and FRMS oc T^ 1 / 2 . For each cases, both 
FRMS and FRMS (Eqn. 29) took their minimum values at r = 0.014 s. The 
resolution of the grid on which I computed RMS prediction errors is 0.001 s 
in the neighborhood of the 0.014 s. 

In a second simulation, I set T = 400 s and A = 100 Hz. I vary r so that the 
expected number of trapped protons per run, Ar, varies from .001 to 4. For 
these cases, FBI AS and FRMS closely track the actual values of FBI AS 
and FRMS (Table 2). For the smallest values of r, the accuracies of the 
approximations are highest. I attribute the slight degradation of approxima- 
tion accuracy at the largest values of Ar to the fact that a w (Eq. 31) is a 
monotonically increasing function of r. 

For convenience, I express the fractional RMS prediction error of A as 



where T is the total time of the experiment. The parameter T* is a function of 
A, r and 5. In proton trap neutron lifetime experiments, A depends on exper- 




(32) 
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imental details including the length of the trap; trapping efficiency; and the 
neutron flux (Refs. 2-5) For the cases considered here, I compute T* directly- 
using Eqns. 9-14 and Eq. 32. For fixed values of A, 5, and r, the derived value 
of T* is approximately the same for all values of T. An exception to this rule 
is when S is large and there are very few bins. I attribute this to truncation 
effects associated with rounding N run to an integer (Eq. 1). Thus, one can 
compute T* from simulation data corresponding to one sufficiently large value 
of T and predict FRMS at other large values of T. As a caveat, for very short 
experiments, the asymptotic theory may not apply and a direct simulation 
may be necessary. 

For the case where the dead time 5 is fixed, T* varies as a function of both A 
and r (Figure 2) in a complicated manner. To clarify resutls, I scale T* and 
t by the true trapping rate A (Figure 3). I define T opt to be the value of r 
that minimizes FRMS. Based on Figure 3, the most elucidating way to find 
the optimal data collection strategy is to minimize AT* as a function of Ar. 
For the cases shown in Figures 2 and 3, I conclude that \T opt increases as A 
increases. In a second simulation experiment, I consider cases where A is fixed 
but the dead time 5 varies from case to case. For these cases, as 5 increases, 
so too does Xr opt (Figure 4). 

5 Discussion 

Earlier I stated that the subsample restriction has no practical effect on data 
collection for neutron lifetime experiments of current interest. To make this 
claim more concrete, I compute the probability of observing n + = N run in the 
full sample for the cases listed in Table 1. For these cases, r = 0.014 s and A = 
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1 s 1 and P(N run ) « 10 L8569Ar ™™. Hence for an experiment of total duration 
100 s, P(N run ) « 10- 13169 . 

In the study, I quantified FBI AS given knowledge of A and particular values 
of r, 5 and T. In actual experiments, one would use the estimated value of A 
rather than the true value. Hence, in Eqn. 28, one would use p = exp(— Ar) 
rather than the true value of p . If FBI AS is negligible, there is no need 
to correct A for bias. In principle, when bias is significant, a bias-corrected 
maximum likelihood estimate should be obtained using the following iterative 
procedure: 

\k+i,BC) = — (33) 

l + FBIAS(\ k , BCh N run ,r) 

where A( ,.bc) — A and \k+i,BC) is the bias-corrected maximum likelihood 
estimate at the fcth iteration. In practice, one iteration of the above proce- 
dure may yield a numerically stable estimate of the bias-corrected maximum 
likelihood estimate for cases of interest. 



6 Summary 

In this work, I studied the statistical properties of a maximum likelihood 
estimate of the rate at which protons are trapped. This study is relevant to 
proton trap neutron lifetime experiments at NIST and similar experiments 
elsewhere. After the first proton is detected, the detector goes dead. Hence, 
the detector provides incomplete information. Due to this incompleteness, I 
showed that the first two moments of the maximum likelihood estimate of the 
trapping rate A are infinite. Hence, the variance of the maximum likelihood 
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estimate is not denned. To construct a maximum likelihood estimate with a 
finite variance, I restricted attention to a subsample of realizations of the data 
that excludes an exceedingly rare realization of the data that yields an infinite 
valued estimate of A. I demonstrated that the probability of observing this 
rare realization quickly decreases to a negligible value for a judicious choice 
of the trapping time for proton trapping rates achievable at NIST (Secion 
5). Hence, restricting attention to the subsample of interest has no practical 
effect on current neutron lifetime experiments of interest. Based on the discrete 
probability density function for this subsample, I derived exact formulas for the 
first two moments of the maximum likelihood estimate of A (Eqns. 10 and 11). 
I derived asymptotically valid formulas for the fractional bias, fractional RMS 
prediction error and fractional standard deviation of the maximum likelihood 
estimate (Eqns. 28-30). I showed that the approximation error associated with 
these formulas is low for a variety of cases (Tables 1,2). 

I demonstrated that the fractional bias (FBI AS) of the estimate was more 
sensitive to r than the fractional mean-square prediction error (FRMS) was 
(Figure 1). As as a function of total observing time T, I showed that FBI AS 
decreases much faster than does FRMS (Table 1). 

I presented an objective method to select the optimal value of r by minimizing 
FRMS. In general, the optimal trapping time r that minimizes FRMS is a 
complicated function of both dead time 5 and the trapping rate A (Figures 2- 
4). For experimental planning purposes, my asymptotic approximations (Eqns. 
28-30) should be useful for determining the optimal data collection strategy 
and for quantifying random and systematic errors. 

In this study, I neglected physical sources of systematic error due to effects such 
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as a time varying proton trapping rate or fluctuations in the actual trapping 
stage interval about the nominal value sought by the experimenter. Hence, 
the bias I quantified here is a purely statistical artifact due to the fact the 
maximum likelihood estimate of the trapping rate is a nonlinear function of 
the observed data. 
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Table 1. Proton trapping rate is A = 1 s _1 . Dead time is 5 = 0.0001 s. Trapping 
stage duration is r = 0.014 s. FRMS and FBI AS are fractional bias and fractional 
root-mean-square prediction error of the maximum likelihood estimate of A. The 
approximations for FRMS and FBI AS are computed using Eqns. 28 and 29. 



T(s) 


N run 


E(n+) 


FBIAS 


FBIAS 


FRMS 


FRMS 


10 


709 


9.85684 


0.000710859 


0.000710858 


0.318749 


0.318742 


25 


1773 


24.6491 


0.0002841 


0.0002841 


0.201479 


0.201477 


50 


3546 


49.2981 


0.000142023 


0.000142023 


0.142446 


0.142446 


100 


7092 


98.5962 


7.10046e-05 


7.10046e-05 


0.100717 


0.100717 


200 


14184 


197.192 


3.55006e-05 


3.55006e-05 


0.0712154 


0.0712154 


400 


28368 


394.385 


1.77499e-05 


1.77499e-05 


0.050356 


0.050356 
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Table 2. Simulation study. Proton trapping rate is A = 100 s . Dead time is 5 = 
0.0001 s. T = 400 s. 



r (s) Ar N run E(n+) FBI AS FBI AS FRMS FRMS 

le-05 0.001 3636360 3634.55 1.37569e-07 1.37569e-07 0.0165873 0.0165873 

le-04 0.01 2e+06 19900.3 2.51254e-07 2.51254e-07 0.00708878 0.00708878 

0.001 0.1 363636 34604.5 1.4461e-06 1.4461e-06 0.00537793 0.00537793 

0.01 1 39603 25033.9 2.1695e-05 2.1695e-05 0.00658726 0.00658698 

0.015 1.5 26490 20579.3 4.38173e-05 4.38173e-05 0.007644 0.007643 

0.02 2.0 19900 17206.8 8.02887e-05 8.02887e-05 0.00896179 0.00895891 

0.03 3.0 13289 12627.4 0.000239664 0.000239664 0.0126488 0.0126306 

0.04 4.0 9975 9792.3 0.000674739 0.000674717 0.0184136 0.0183143 
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Total Observing Time = 50 s 
X = 1 s" 1 Dead Time = 0.0001 s 
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Figure 1. Fractional bias (FBI AS) and fractional standard error (FSE) for 
simulation experiment where true trapping rate is A = 1 s _1 , total duration of 
experiment is T = 50 s, dead time is 5 = 100 /US, and trapping stage duration 
r varies. 
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Dead Time = 0.0001 s 
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Figure 2. Fractional RMS prediction error (FRMS) of the trapping rate is 
expressed as O.OOly^r, where the total length of the experiment is T. Dead 
time is 5 = 0.0001 s. 



19 




Figure 3. Same results as in Figure 2, but T* and r are rescaled to clarify 
results. 
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